Energy non-equipartition in systems of inelastic, rough spheres 
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We calculate and verify with simulations the ratio between 
the average translational and rotational energies of systems 
with rough, inelastic particles, either forced or freely cooling. 
The ratio shows non-equipartition of energy. In stationary 
flows, this ratio depends mainly on the particle roughness, but 
in nonstationary flows, such as freely cooling granular media, 
it also depends strongly on the normal dissipation. The ap- 
proach presented here unifies and simplifies different results 
obtained by more elaborate kinetic theories. We observe that 
the boundary induced energy flux plays an important role. 
PACS: 51.10.+y,46.10.+z,05.60.+w,05.40.+j 

Granular materials are collections of macroscopic par- 
ticles with rough surfaces and dissipative interactions, 
as addressed in this letter. Although rotation and fric- 
tion are often neglected, they play an active role for the 
dynamics of systems with rough or non-spherical con- 
stituents. In contrast to classical elastic systems, energy 
is not equipartitioned between the degrees of freedom in 
the system |l|-|5| . In order to examine this ratio, kinetic 
theories [[jj-Q] and numerical simulations |5| were applied 
for special boundary conditions, and a variety of results 
were obtained. We unify these results in a single theory 
which also explains when each one is valid. 

We consider a system of N particles. We define E to 
be the average translational kinetic energy per degree of 
freedom, and E° to be the rotational kinetic energy per 
degree of freedom, so that: 
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Here, m is the mass and a the radius of a particle. Vi is 
the velocity of particle i, and uji is its angular velocity. 
T> is the number of dimensions (we restrict ourselves to 
T> = 2 and 3 here), q is the dimensionless moment of 
inertia; q = 1/2 for disks, and q = 2/5 for spheres. The 
number of translational degrees of freedom per particle is 
n = T>, and the number of rotational degrees of freedoms 
is n° = 2T> — 3. E and E° are often referred to as "gran- 
ular temperatures" . This terminology is not intended 
to suggest that a thermodynamic equilibrium exists in 
granular flows, but simply to draw an analogy with the 
temperature of an ideal gas, which is also the average 
energy per degree of freedom. E and E° in Eq. (|l|) are 
well defined whether or not the system is in equilibrium. 
In this paper, we consider all particles to be identical, 
however, the above definitions can easily be extended to 
different types of particles. 

We quantify the distribution of energy between the 
translational and rotational modes with the quantity 



R = E° I (E° + E). When there is no energy in the rota- 
tional mode, R = 0. When energy is equally distributed 
between all the modes, R = 1/2. If rotational energy 
dominates, then R — > 1. We will study how R depends 
on the particle properties and the boundary conditions. 
This question has been addressed by several authors |l]-|4|] 
but serious and unexplained conflicts exist between their 
results. For example, Q claims that R depends on the 
normal restitution, whereas 0H^| say R is independent of 
this parameter. 

We use the standard constant roughness model for the 
instantaneous collisions of rotating particles with radius 
a, mass m, and moment of inertia I — qma 2 . This model 
accounts for dissipation, using the restitution coefficient 
r and the tangential restitution j3. Since it has been 
extensively used and discussed Jl|-|| , we include only the 
results here. The post-collisional velocities v', uj' are 
given in terms of the pre-collisional velocities v, u> by 
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Here, v n = [(vi — V2) ■ n] ■ n is the component of v-y — 
parallel to h, a unit vector pointing along the line con- 
necting the centers of the colliding particles. The tan- 
gential component of t>i — V2 is v t = f 1 — V2 — v n and 
v r = — a(u>i + o>2) x n is the tangential velocity due to 
particle rotation. 

Later on, we will need expressions for the change in 
rotational and translational kinetic energy during a col- 
lision. The change in translational energy is 

NnAE = -Qv 2 n + S [~C tl v 2 t - C t2 (v t ■ v r ) + C t3 v 2 ] , (3) 



with the positive prefactors Q = m(l—r 
/3)/[4(l + q) 2 ], and the constants C t \ 
C t 2 = 2 — 2q(3 and C t z = q(l + (3). Likewise, the change 
in rotational energy is 



2 )/4, S ee mq(l + 
= 2 + g (l-/3), 



Nn°AE° ee +S [C rl v 2 - C r2 (v t ■ v r ) - C r3 v 2 ] 



(4) 



where the constants are C r \ = (1 + /?), CV2 = 2(q — (3), 
and CV3 ee 2q + 1 — j3. Note that the C are positive 
(except that C r i can be negative) so that the signs in 
Eqs. (||) and (^) indicate the direction of energy transfer 
between the degrees of freedom. 

The quantity R (or its equivalent) has been calculated 
by several authors. One result found by three different 
authors for granular material undergoing uniform 

shear is 
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R = 



q(l + (3) + (1 + 2q - P) ' 



(5) 



i.e. a function independent of r. 

On the other hand, Goldshtein and Shapiro Q found 
a much different expression (for T> = 3), involving 7': 
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with the quantities ac = (1 — P 2 ){\ — + q) — l + r 2 , 
and b G = 2q{\ + /3) 2 /(l + q) 2 . Eq. (§) differs greatly from 
Eq. (|J) when r < 1. 

In the following, we show that the differences between 
Eqs. (§) and (§) arise from the boundary conditions, i.e. 
from the existence of external forcing. 

Consider a granular material with external forcing, 
where the particles interact only through collisions obey- 
ing Eq. (||). The change in E° from time to to t\ is 

E°( tl ) - E°(t ) = AE °&) + P °(t) 



coll. 



(7) 



where P° is the rotational energy added by the forcing. 
The sum is taken over all the collisions which take place 
between to and tx, and AE°(Ci) is the change in E° for 
collision i. Now, consider a situation where the granular 
medium is maintained in a stationary state [E°(ti) — 
E°(to)} by some kind of forcing, and that this forcing 
adds only translational kinetic energy [P° — 0]. Eq. (Q) 
becomes: J2 C oii AE°(d) = 0, which states that collisions 
do not, on average, change the rotational energy. When 
the assumptions made above are satisfied, the dissipation 
of rotational energy is exactly balanced by the conversion 
of E into E°. Using Eq. (@) we get: -C* r i(v 2 ) + C r2 (v t ■ 
v r ) +C r z{v 2 ) — 0. Here, the angle brackets (...) indicate 
an average taken over the collisions. We now consider 
how to calculate these averages. 

If if; is the quantity to be averaged, then 



(tp) = {i>P co ]\(Vl,V2,0>l,U2,b)}b, 



(8) 



where P C oii gives the probability of a collision occuring 
between particles 1 and 2 with a normalized impact pa- 
rameter < b < 1. The normalized impact parame- 
ter is the distance between particle centers at closest 
approach if the particles did not interact, normalized 
by the particle diameter. The subscripts on the brack- 
ets means we average over all values of b, and over all 
pairs of particles. We now make several simplifying as- 
sumptions about Pcoii- First, we assume that the an- 
gular velocities have no effect on the probability of col- 
lision: P co ii = P co i\(v l) i>2, b). Next, we assume that 
Poll ~ vf(b), (v = \vi — v 2 \) because particles with large 
relative velocities are more likely to collide. Finally, the 
dependence of P co ii on b can be deduced from geometrical 
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with the factor of 2 and the denominator required for 
normalization. 

Since the rotational and translational velocities are un- 
corrected, 



?) = {^} 1 , 2 = 2 a 2 {(u, 1 xn) 2 } 1 = 
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The factor of T> — 1 arises because in 2? = 3, one of 
the three rotational degrees of freedom is excluded by 
the cross product with n. In T> — 2, there is only one 
rotational degree of freedom, and it always participates 
in every collision. 

To calculate (v 2 ) and (vf), we use vf = v 2 b 2 and v 2 = 
v 2 — v 1 . The average can be factored into two parts: 
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Evaluating the first factor gives {b 2 (2b)^- 2 } b = D/6. 
The second factor will be proportional to E/m, but calcu- 
lating the coefficient requires knowledge of the distribu- 
tion of velocities. Assuming a Maxwellian velocity distri- 
bution gives {v 3 } h2 /{v} h 2 = 24(2? - l)E/(Vm). Then, 
we have (v 2 ) = 4(P - l)E/m and (tj 2 ) = 8E/m. 

All of the assumptions we have made up to know are 
equivalent to those made in the kinetic theories 
Thus, it is not surprising that we recover some of their 
results. In particular, putting the averages into Eq. (0) 
gives ~C r \E + C r zE° jq = 0, and after using the defini- 
tions of C r \ and CV3, we obtain R as in Eq. (||). 

Eq. (||) does not depend on r because R is fixed by 
a balance between the conversion of E into E° and the 
dissipation of E°. Both of these processes depend only 
on j3 and q but not on r. As soon as the dissipation of E 
starts to play a role in determining P, then r will appear. 
We examine such a case next. 

Consider a granular medium in the absence of forcing. 
If r 7^ 1 and j3 7^ 1, — 1, then E° and E decrease towards 
0, but R can approach a constant. This can be verified 
by simulations and a more elaborate calculation in the 
framework of the kinetic theory Q or with a Liouville 
operator formalism ||. 

In the following, we simplify the algebra by using K — 
E/E° instead of P [P = 1/(1 + K)]. During a collision, 
K changes by 



AK 



E + AE 
E° + AE° 



K 



AE - KAE° 
E° + AE° 
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We look for a value of K such that AK = 0. The de- 
nominator of this equation is always positive. Equat- 
ing the numerator to 0, we find that a collision leaves 
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R unchanged if K = AE/AE°. This equation can be 
expanded in terms of (v„), (v 2 ), and (v%) using Eqs. (||) 
and (|), to 

A(vl)+C tl (v?)-C t3 {vi) 
-aC rl (v 2 } + aC rS (v?) ' 1 ' 

where a = n/n°, and A = Q/S [see Eq. (j|)]. Because the 
energy decreases with every collision, the averages must 
be interpreted as taken over all possible collisions at a 
given time. 

Using our previous expressions for the averages and 
reorganizing Eq. (|l3|) as a quadratic equation for K yields 

a-^-K 2 -[aa G + (a-l)c G ]K - y = 0, (14) 

with the quantities a G and b G from Eq. (^), and c G = 
(1+ (3)(2q + q 2 - (3q 2 )/(l + q) 2 . In deriving Eq. (0), we 
used a = 2/(23 — 1). For T> = 3, we have a = 1, and the 
solution of Eq. (Q) leads to Eq. (§). 

Next we compare the theoretical results, derived above, 
with simulations in D = 2. We examine three different 
simulational "experiments" . In the first case, energy is 
put into one translational mode by a vibrating wall, and 
in the second case, a granular material under sheared pe- 
riodic boundary conditions is examined. Finally, a gran- 
ular media is studied in the absence of any forcing what- 
soever. On the basis of the theory presented above, we 
expect the first two cases to obey Eq. (|5|), and the last 
case to obey Eq. (^]). In all cases, we perform the exper- 
iments with different r, and for each value of r, we vary 
/3 from -0.95 to 1. 

In the first experiment N — 160 particles of radius a 
are placed on a vertically vibrating floor in the presence 
of gravity. The boundaries in the horizontal direction 
are periodic, the domain is 50 particle radii wide and in- 
finitely high. The period of the floor vibration T and the 
gravitational acceleration are here related by gT 2 /a = 1. 
The height of the floor varies periodically in time, follow- 
ing an asymmetric sawtooth wave form. (The choice of 
wave form is arbitrary; changing the wave form does not 
significantly change R.) The floor moves up a distance 
of 5a, with upwards velocity 5a/T and then returns in- 
stantly to its lowest position. In all cases, the simulations 
ran for 7000T, with R being measured every At = 0.25T 
for 5000T < t < 7000T, and these values were aver- 
aged to give the points in Fig. 0(a). Since this experi- 
ment satisfies the assumption that energy is input into 
the translational modes only, we expect that the results 
will satisfy Eq. (||). Fig. 0(a) confirms that this is indeed 
the case, besides some systematic underestimation of R 
by the theory when the dissipation is strong. 




-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 



p p 

FIG. 1. Simulation results from a vertically vibrated sys- 
tem, (a) R is compared with the theoretical prediction in 
Eq. (0) for three values of r. (b) The data for r — 0.9 from 
Fig. nfa) are compared to simulations with different boundary 
conditions (see text for details). 

Since it is difficult to do experiments with periodic 
boundaries, we also examined the effect of stationary side 
walls at the edges of the domain. We compare the re- 
sults for r = 0.9 from Fig. pl(a) with three different types 
of boundary conditions. (!) a perfectly rough wall with 
r — [3 — 1, (ii) a wall where r and /3 have the same val- 
ues for particle-particle and particle-wall collisions, and 
(iii) a perfectly elastic and smooth wall with r = 1 and 
[3 = —1. The results presented in Fig. 0(b) show that 
only the first type of wall causes R to deviate signifi- 
cantly. We relate this to a competition between particle- 
particle and particle-wall collisions: The collisions with 
the wall in case (i) push R towards R{(3 = 1) = 1/2, 
i.e. equipartition holds, and the collision between par- 
ticles push R towards a smaller value. (At (3 P = — 1, 
R > 1/2 because the kinetic energy of the vertical veloc- 
ities is greater than that of the horizontal, and the wall 
couples only to the vertical motions.) In case (ii), both 
types of collisions push R towards the same value, and 
in case (iii), the particle- wall collisions do not influence 
R, so that the results are not perturbed at all. 

In the second experiment, we drive the granular ma- 
terial by shearing it, a case frequently considered in the 
literature |l]-||]. N — 160 particles are placed in a square 
domain whose sides are L — 50a in length. The bound- 
aries are periodic in the x and y direction. A uniform 
shear is imposed by applying Lees-Edwards boundary 
conditions H: when a particle exits the domain at the 
bottom (top), its image enters at the top (bottom) of 
the domain, with its x velocity increased by a constant 
velocity U (-U), at its position shifted by a distance Ut 
{—Ut). These boundary conditions eliminate the need to 
specify wall properties and make the system translation- 
ally invariant. In our simulations, U ~ 2a/T, giving a 
shear rate of T = U/L = 0.04/T. The time unit T is 
arbitrary. 

A direct application of Eq. (|B|) will fail, because the 
boundary conditions generate both an average flow and 
an average rotation. However, if we interpret E° and 



3 



E to be the energy which remains after removing the 
mean flow and rotation, the agreement between theory 
and simulation is good, as shown in Fig. §(a). The trans- 
lational kinetic energy of the mean flow was calculated as 



-Bmean = (m/2)r 2 E^ife" L I 2 ? The rotational energy 
of the mean flow was estimated by fi^,, = mqa 2 rt 2 /2, 
where f2 is the observed average angular velocity. Even 
with longer averaging times (5000T in our simulations), 
the data is considerable noisier than in Fig. [|. 



eludes boundary conditions with with non-zero rotational 
energy input, and calculating R for the more realistic in- 
teraction model that accounts also for Coulomb-friction. 
Ill 

We acknowledge the support of the "Alexander von 
Humboldt-Stiftung" and of the "Deutsche Forschungsge- 
meinschaft, Sonderforschungsbereich 382" . 



0.6 



0.4 



0.2 



0.0 



Eq. (5) 




* r=1 .0 




o r=0.9 




o r 0.8 




□ r=0.7 j 






(a) 



-1.0 -0.5 0.0 



0.5 




FIG. 2. (a) Results of the shearing experiment. The 
points show the results for different r, compared to the an- 
alytical result, (b) Results of the cooling experiment. The 
points show the results of simulations and the curves give 
R = 1/(1 + K) with K from Eq. (0), with the appropriate 
value of r, and D — 2. 
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In the last experiment, N = 160 particles are placed 
in a periodic domain with L = 100a, but no shear is ap- 
plied. The initial condition is generated by setting r = 1 
and j3 = — 1. The particles quickly attain a Maxwellian 
velocity distribution. Then, dissipation is "switched on" , 
and the system evolves without any further input of en- 
ergy. Although the energy decreases with every collision, 
R approaches the constant value shown in the graph. The 
results of many simulations were averaged together to re- 
duce the fluctuations. In Fig. ^(b), we plot the results. 
The theoretical curves are the solutions of Eq. (|i~i|), with 
a = 2, and the three values of r shown in the plot. 

In conclusion, we extend Eq. (^J) to two dimensions 
and summarize the different results for R in the litera- 
ture. Our method of calculation is simple enough to show 
why Eq. (^|) applies to forced granular media, and Eq. (|^) 
to cooling granular media. Eq. (|5|) will apply whenever 
the forcing adds only translational energy. This is the 
case for both vibration and shear. We succeeded to for- 
mulate a procedure to calculate the ratio of tangential 
and rotational energy in dissipative systems. Our result 
replaces the value K = 1/2 (in 2D when equipartition 
is true for elastic systems). We examined not only the 
limit of almost elastic particles, but also rather inelas- 
tic situations. For most boundary conditions used, the 
agreement between theory and simulations is encourag- 
ing. However, simulations in T> = 3 are stil needed to 
check the analytical expressions. Possible future work in- 
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